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ABSTRACT 

We investigate the launching of outflows in the close vicinity of a young stellar object, treating the 
innermost portion of an accretion disk as a gravitationally bound reservoir of matter. By solving 
the resistive MHD equations with our version of the Zeus-3D code with implemented resistivity, we 
study the effect of magnetic diffusivity in the magnetospheric accretion-ejection mechanism. Physical 
resistivity has been included in the whole computational region. We show, for the first time, that 
quasi-stationary outflows consisting of axial and conical components can be launched from a purely 
resistive magnetosphere. We identify four stages of magnetospheric interaction with distinctly differ- 
ent geometries of the magnetic field, and describe the effect of magnetic reconnection in re-shaping the 
magnetic field. The stages are the relaxation, reconnection and infall, after which two outflow compo- 
nents can be seen in a final flow: a fast axial component launched from above the star, dominated by 
magnetic pressure, and a slow conical component, launched from the opened resistive magnetosphere 
of a disk gap, between the star and the disk inner radius. We show how outflows depend on the 
disk to corona density ratio and on strength of the magnetic field, and compare the position of the 
disk truncation radius with theoretical predictions. Results from previous investigations with resistive 
MHD in the literature, which have been obtained with various setups, are recovered in our simulations. 
Comparisons are thus made easier for more general purposes, by identifying previous features in the 
simulations within the different stages of our simulation. 

Subject headings: methods: numerical — processes: MHD — stars: formation 



1. INTRODUCTION 

Highly collimated outflows have been observed from 
AGNs t o Young Stellar Obje c ts (YS Os) and young brown 
dwarfs (jWhelan et all 120051 120071) . Accreting compact 
stars, like accreting white dwarfs in symbiotic binaries 
( Sokoloski fc Kenvonll2003ft and neutron stars like Cir X- 
1 (jHeinz et al.1 12007ft . also show similar outflowing phe- 
nomena. An outflow is characterized as a jet if it is super- 
magnetosonic, collimated into an apparent narrow open- 
ing, and reaches a stationary or quasi-stationary state. 
Such high- velocity outflowing fluxes of matter are an in- 
tegral part of stellar evolution. Observations in multiple 
wavelengths are reaching closer and closer to the objects 
that drive them. 

Among all the systems, models of launching outflows 
in YSOs are closest to scrutiny by observations due to 
available data from star-forming regions. An accre- 
tion disk, through which matter accretes onto the young 
star with velocities close to a free-fall, is often associ- 
ated with a jet-driving YSO (|Edwards et al.lll994 1^555 ) . 
A correlation between the accretion rate and the high- 
velocity jet pow er was found in ma ny Classical T-Tauri 
Stars (CTTSs) (jCabrit et al.lll990f ). The ratio of mass 
loss in the outflow to disk accretion rate, M w /M a , ex- 
tracted from observations is hard to constrain. It is 
best estimated to be approximately 0.1 through mea- 
sure ments of optica l forbi dden l ines and v eiling — see 
e.g. lHartigan etHI (|1995f ) and lEdwardsl (|2008l ). Re- 
cently, He I A 10830 has offered a good probe into the 
high-velocity winds originating from the inner region 
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where the star interact s with the disk (|Edwards et al.1 
I2003t iKwan et al.|[2007h . Despite the potential diagnos- 
tic power of such emission lines, the actual structure and 
physical conditions of outflows can be more complex. To 
further interpret the observed line profiles, and the ori- 
gins of outflows from the close vicinity of CTTSs, pre- 
dictions from both theoretical and numerical models are 
required. 

Outflows driven by energy derived from accretion are 
particularly appealing in the scenarios of jet launch- 
ing. Many models have been proposed based on 
the concept of magnetoc entrifugal wind mechanisms 
(jBlandford fc Pavnd 11982) . differing in the origins of 
the underlying magnetic fields and locations of matter 
launching. An outflow could be a disk wind driven by 
magnetic fields dragged in from the envelope or gener- 
ated by the disk dynamo, or an inner disk wind an- 
chored to the narrow innermost region as in the X- 
wind model powered by an enhanced dynamo from 
the star-disk interaction (|Shu et all 11994 11997ft. si- 
multa neously with an accretion funnel (|Ostriker fc Shul 
1995). It might also be a stellar wind driven along the 
open field lines fro m the stellar surface by thermal or 
magnetic pressure (Ivqn Rekowski fc Brandenburg! 120061 : 
iRomanova et al.l 12005ft . or some combination of the dif- 
ferent possibilities. Related to the launch of winds, mag- 
netospheric accretion has been des cribe d in works by 
Konig] dl99l . lOstriker fc Shul (|l995ft and lKoldoba"eTaLl 
(2002) in the context of a magnetosphere interacting 
with the surrounding disk, sharing some similarities with 
the co mpact objects like neutron stars (|Ghosh fc Lamb! 
Il979afbh . Except for the pure disk wind models, a mag- 
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TABLE 1 

List of assumptions for initial conditions in some relevant works 



Paper 








LUi Wllct 


Havashi et al. (1996) 


10 


non-rotating 


in rotational equilibrium & adiabatic 


isothermal, non-rotating 


riirosc ct al. (iyy/j 


1U 


non-rotating 


adiabatic, Keplerian 


isothermal, hydrostatic 










rotates different than disk 


Mil lor JCj- Cltnnp M QQ7^ 


i 


rotating 


adiaDatic, Ivcplcrian 


isothermal, solid body 










corotating with star at R COT 


Romanova et al. (2002) 


10 2 


rotating 


adiabatic, super-Keplerian 


adiabatic, corotating with star 










for R< i?cor, else with disk 


Kiiker et al. (2003) 


10 4 


rotating 


adiabatic, Keplerian 


not in hydrostatic balance, 










non-rotating 


Ustvuerova et al. (2006) 


10 3 


rotating 


adiabatic, sub-Keplerian 


adiabatic, corotating with star 








for R< Rcor, else with disk 


Romanova et al. (2009) 


10 4 


rotating 


isothermal, sub-Keplerian 


isothermal, corotating with star 










for R< i?cor, else with disk 



netically connected star-disk system plays an important 
role in the making of the young stellar system and the 
evolution of angular momentum through the generation 
of strong outflows during the main phase of accretion. 

Numerical investigations have been followed up on the 
time-dependent evolution of a system where the central 
star is magnetically connected to its accretion disk and 
their connection to jet for mation and acc r etion. In one 
of the earliest attempts by lHavashi et al.l (|1996fl , where 
simulations of only a few rotation periods were obtained, 
a dipole magnetosphere corotating with the central star 
threaded the accretion disk that was in Keplerian rota- 
tion. Magnetic field lines connecting both the disk and 
the star inflate outwards due to shear, and reconnection 
blows out the matter along with the field, partially open- 
ing up the originally closed dipole loops. Gas can out- 
flow from those opened field lines and might form part 
of the X-ray jet that is often associated with flares. Re- 
connection as a possible origin o f X-rays from such sys - 
tems has also bee n indicated in Idal Pino et al.l (|2010( L 
iHirose et al.l ()1997f ) investigated a magnetized star inter- 
acting with a truncated disk that was threaded with an 
initially uniform field dragged in from the outer core, 
in the same direction of the magnetosphere, but sepa- 
rated by a neutral current sheet in the equatorial plane 
as a result of interaction between the fields brought to- 
gether. For simplicity, the star was not rotating, but the 
differentially rotating disk could anyway provide enough 
shear to make the field inflate outwards, followed by a 
reconnection event and mass transfer onto the magne- 
tosphere. The transferred mass diverted into two direc- 
tions: one that falls onto the star and the other that flows 
out alon g the opened st el lar fie ld lines. Longer simula- 
tions by iGoodson et all (|1997l ) with an aligned dipole 
and a conducting accretion disk showed that differential 
rotation of the disk can drive episodes of loop expan- 
sion. Such expansion can drive two outflow components 
of gas: one hot convergent flow along the rotation axis, 
and another, slower cold flow on the disk side of the ex- 
panding loop. Mi ller fc Stone! ([T997), on the other hand, 
investigated interactions of magnetospheres with accre- 
tion disks under three different magnetic configurations 
and their respective dynamical evolution. All of the men- 
tioned works involve resistive MHD and simple models 
of accretion disks. 

Numerical treatments of the physical processes and 



disk structures h ave been improved over the years for 
star-disk systems. iKuker et al.l (|2003t ) solv ed the disk i n 
ID with a radiative hydrodynamic code by iKlevI |l989), 
and then extrapolated the solution to 2D as their initial 
condition. For the full 2D axisymmetric MHD problem, 
the induction equation, Lorentz force and Ohmic dis- 
sipation were now included into Kley's code, with the 
assumption of equal viscous and resistive dissipations. 
The main result was that for a smaller magnetic field 
than 1 kG the disk is not disrupted; but for a larger field 
of the order of 1-10 kG, an outflow could be launched 
from the disk. Without reaching a steady state, the cen- 
tral star was spun up by the prevalent angular momen- 
tum transported to it, while the magne tic fi eld acted to 
slow i t down. In lRomanova et all (|2002f ) and lLong et al.l 
( 2005), a star and part of the magnetosphere corotated 
up to the corotation radius, and the magnetosphere coro- 
tated with the disk farther out. The corona was treated 
in the ideal MHD regime, with effective numerical re- 
sistivity diffusing magnetic field in the radial direction. 
They found funnel flows onto the central object, spinning 
up or down the star, depending on the ratio of rotation 
r ate of the star to the rot ation rate of the disk inner rim. 

Rom anova et al.l (|2009l ) (hereafter R09) later investi- 
gated the effects of physical viscosity and resistivity. 
When the magnetic Prandtl number, the ratio between 
the viscosity and the resistivity Pr = v /rj, is greater than 
one, viscosity can strongly influence the solution. They 
found that, in addition to the fast and light axial jet 
above the star, there is another, new conical wind flow- 
ing up to 30 percent of the matter from the innermost 
portion of the disk. However, they required two different 
setups in their simulations, one for a slow and the other 
for a fast rotating star. In one of them, simulations were 
started with a slowly rotating star without any matter 
in the computational box, and then the stellar rotation 
was gradually speeded up to its maximum value, with 
matter slowly inflowing from the outer boundary. Ini- 
tial relaxation of the interaction between magnetic field 
and matter in such setup was different from most other 
simulations. The comparison with previous results was 
then complicated even more when different initial and 
boundary conditions were used. Also, previous results, 
which did not show a stationary conical outflow in the 
literature, were in the regime of Pr < 1. It was not clear 
if such a component emerged only for a slowly rotating 
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star with viscosity larger than resistivity. 

We aim here to put previously published results in per- 
spective, in systems of magnetospheric star-disk interac- 
tion. Table [1] lists the kinematic and thermodynamic 
assumptions in published works, each of them being usu- 
ally repeated with a variety of parameters or methods. 
Here, we focus on the hydro-magnetic part of the mecha- 
nism, assuming that thermal and radiation pressure bal- 
ance each other in the innermost star-disk region. One 
important parameter to distinguish the models is the al- 
ready mentioned magnetic Prandtl number. In R09 it is 
claimed that if viscosity is larger than resistivity, a new, 
quasi-stationary conical component appears. However, 
we found this component even in our rendering of previ- 
ously studied cases, with resistivity larger than viscosity. 
Another parameter we study is the density ratio between 
the disk and the corona. It is usually included as a free 
parameter of the order of 10 2 or 10 3 , at best 10 5 , without 
much discussion, but from astronomical observations we 
know that this ratio is a few orders of magnitude larger, 
up to 10 8 . We investigate the influence of this ratio on 
the mass and angular momentum flux in the launch of 
outflows. There are other possibilities in the setup, which 
we did not investigate here, e.g. inclusion of the stellar 
wind, which would probably affect the open stellar field. 

Most of the previously published findings in magneto- 
spheric interactions could be identified in our simulations 
within four evolutionary stages of a single simulation, 
with both the star and the disk set in the computational 
box from the very beginning. The magnetic field has 
similar topology in the corresponding stages of different 
simulations, and the final stage differs only in the rela- 
tive intensity of the two outflow components. This in- 
tensity can be influenced by the dissipative mechanisms. 
Resistivity, which controls the onset of magnetic recon- 
nection, triggers the necessary change in the geometry of 
the magnetic field needed for the outflow launching. 

The organization of the paper is as follows. We first 
describe our implementation of the boundary and initial 
conditions. In £[3] we report regimes we found under a 
broad range of parameters. We investigated the influence 
of corona to disk density ratio, strength of magnetic field 
and the physical resistivity. In 2] we address the role of 
reconnection in the launching, in [J5]we check a criterion 
for the site of launching, and in [J6]we compare position of 
the disk truncation radius in our simulations with some 
theoretical predictions. Then we discuss investigated pa- 
rameters and the resulting outflows. 



NUMERICAL SETUP FOR THE RESISTIVE MHD 
SYSTEM 



We extend previous work of ICemeliic feFm dt (2004), 
who adopted a disk in the resistive MHD regime 
and its halo in the ideal-MHD regime, following 
iCasse fc Keppensl (|2002t ). We implement an absorbing, 
rotating stellar surface layer enclosing the origin, and 
include resistivity in the whole computational box. The 
initial conditions of density and magnetic field are shown 
in Figure [TJ and the setup of the stellar surface as a 
boundary layer inside the computational box is shown in 
Figure [H 

The equations of resistive MHD are solved using our 
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Fig. 1. — Standard initial and boundary conditions in our sim- 
ulations. The initial hydrostatic density distribution in the disk 
corona and the disk is plotted in logarithmic color grading. The 
density in the disk is four orders of magnitude larger than in the 
corona. The dipole stellar magnetic field is plotted in white solid 
lines, and velocity vectors are shown in white arrows. The stellar 
surface is defined as a rotating absorbing boundary layer, enclos- 
ing the origin — see the zoom into this region in Figure [2] In 
simulations S2, the stellar absorbing layer extends into the small 
portion of the disk mid-plane inside the disk gap, of radius R g , as 
an outflow boundary. Along the axis of symmetry and at the mid- 
plane in the disk, a reflection and anti-reflection boundaries are 
imposed. An outflow boundary is imposed on the outer bound- 
aries of the computational domain, except for the disk outer rim, 
where a small inflow into the disk is set, to mimic the accretion 
flow from the portion of the disk beyond the computational box. 




Fig. 2. — Zoom into the setup of the stellar surface from the Fig- 
ure [l] The star is set as a rotating, absorbing boundary condition 
inside the computational box, enclosing the origin. Components 
of the poloidal velocity tip are copied from the layer immediately 
above the star. The stellar rotation rate, determined by the initial 
toroidal component of the velocity vt at the stellar surface, is kept 
constant throughout the simulation. 



version of the Zeus-3D cod<E Zeus347 (|Fendt k. Cemeliid 
2002), in axisymmetry option. They are, in the cgs sys- 
tem of units: 

dp 

at 



+ V-(pv) = (I) 



dv 
~dt 



+ (v-V)v 



jxB 



(2) 



1 For general description and n umerical methods used in Zeus 
code sec Stone & Norman ( 1992a, b) 
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dB 



- V x 



47T 

v x B rtf 



= (3) 



<2i 



+ (v • V) e 



•p(V-v)=0 (4) 



3 = t-V x B 

4tt 



(5) 



where we neglected the Ohmic term in the energy equa- 
tion. For a complete set of equations, an ideal gas law 
is assumed. The symbols in the equations of continuity 
of mass, momentum and induction equation have their 
usual meaning: p and p are the matter density and ther- 
mal pressure, v is the velocity, $ = -GM*/(R 2 + Z 2 ) 1 / 2 
is the gravitational potential of the central object, and 
B and j are the magnetic field and the electrical current, 
respectively. In cgs units, magnetic diffusivity is equal to 
resistivity, so r/ stands for the electrical resistivity. The 
resistive term in the induction equation (Equation [3]) is 
included in the code by subtracting Airqj jc from the elec- 
tromotive forces in the MOCEMFS procedure in Zeus347 
— see Appendix A in lFendt fc Cemefiicl (|2002D for tests. 
In the energy equation, e = p/(7 — 1) is the internal 
energy per unit volume. The density is related to the 
pressure by a polytropic and isentropic relation, and the 
initial entropy follows S = ln(p/p 7 ), with an adiabatic 
index 7 = 5/3 in the initial conditions. 

We solve these MHD equations in dimensionless form. 
The variables are normalized to their value measured in 
the mid-plane of the disk, at a fiducial radius Rq, which 
we choose to be inside the initial disk gap, Ro — 2.85i?*, 
where i?* is the stellar radius. We include radial dis- 
tances up to 20 stellar radii in our computational box. 
The actual dimensional values are determined by sub- 
stituting the mass of the central object M„, and the 
fiducial radius and density, with the mass accretion 
rate for the disk, in units of Mq = -Ro/°o v k,Oj where 
v k,o = {GM^/Rq) 1 / 2 is the Keplerian velocity at Rq. 
The normalized coordinates are R' = R/Ro, Z' = Z/Rq, 
and v' = v/vk,o- The time in the code is measured in 
units of the rotation rate timescale at Rq, to = -Ro/vk.o, 
and the period at R is equal to 2n. The dimensionless 
equation of motion can be written as: 



— + (V • V) v' = J — : - ^7 - W 



Mi p' 5 oP > 



(6) 



with V = i? V, t' = t/t , p' = p/po, B' = B/B and 
$' = — \/ [R 12 + Z' 2 ) 1 / 2 . Primes are omitted in equations 
in the rest of this paper, and quantities are written in 
code units, unless otherwise specified. 
We introduced the free parameters: 



M 



A,0 



47F A)Vk / b o and <5o = PoVk q/po • (7) 



The Alfvenic Mach number, Ma.o, at Rq and Z = 0, 
determines the magnetic field strength. In a typical run, 
Ma — 225 for a dipole field of the order of 100 G, at the 
surface of the star. The kinetic to thermal energy density 
ratio, do, is the square of the gas Mach number, whose 
fiducial value can be estimated from the definition of the 
adiabatic coefficient. For our setup we choose 5q = 100 
or 50. The typical temperatures in the corona and the 

2 The sound speed is c 2 = (dp/dp)\g = 'y^R.T/m, where S denotes 



disk of the YSOs are 10 6 K and 5 x 10 3 K, respectively, so 
that at the inner edge of the disk it is c s , corona : c s .disk = 
0.8 : 0.05, in the units of Vk,o- 

We use a uniform grid of R x Z = (90 x 90) cells, 
in the axisymmetric cylindrical coordinates (R,4>,Z). 
The physical scale corresponds to (20 x 20) stellar radii. 
We performed larger domain simulations in R x Z = 
(60 x 60)i?*, and simulations in higher resolution with 
RxZ = (180 x 180) grid cells in (20x20)i?*. These simu- 
lations lasted long enough for comparisons, but are more 
prone to numerical problems, and tend to cease during 
the relaxation or soon afterwards, so that it is harder to 
perform a thorough parameter study using them. This is 
probably because of numerical viscosity, which is larger 
with less resolution, and helps the code to go through 
problematic events. Inclusion of physical viscosity would 
enable larger resolution, but we focus here only on effects 
of resistivity. 

2.1. Example of rescaling 

We give an example of rescaling for a case of an YSO 
with M* = 0.8M Q , i?* = 2i?0, so that fiducial distance is 
Ro = 5.7-R© = 0.027 AU. Our computational domain is 
then R x Z « 0.2 x 0.2 AU. We can rewrite the Keplerian 
speed at Rq in units of solar mass and radius as 



vk,o 



GM* 



<GMr. 



M* / (Ro_ 



Rq y Mq \R 



(8) 



giving the fiducial velocity 1.64 x 10 7 cm s _1 . The period 
of Keplerian rotation at Ro is Po = 2ttRo/vk,o = 1-76 
days. The stellar rotation rate for T-Tauri Stars is 
usually about 1/10 of the breakup rate, which we ob- 
tain from {GAh/Rlf/ 2 = 2 x lO^s" 1 = 0.4 days. 
This means that period of rotation of a star should 
be about 4 days. If we assume an accretion rate of 
Mq = 10 _8 M Q yr _1 , fiducial density and pressure are 
po = 2.44 x 10~ 13 g cm -3 and p Q = 109 erg cm -3 , re- 
spectively. With c s ~ vk,0i the reference temperature is 
T = mv K o/(7^) ~ 10 6 K- The reference value of re- 
sistivity is 770 ~ vk,o-Ro = 10 19 cm 2 s _1 , which is much 
larger than the classical Spitzer value. 

The strength of the magnetic field we obtain from the 
magnetic pressure at the mid-plane of the disk, Bq(Z = 
0) = (inpoSo) 1 / 2 /Ma,o (see Equation [7]) , and the stellar 
dipole field at Ro is B — S*(i?*/i?o) 3 - A complete 
expression for the fiducial magnetic field we can write in 
units of solar mass and radius as: 



Bi 



47T 



M a M Q GMq M»/M 



M A M Q /yr yr 



R 5 



(Ro/ Rq 



(9) 



The factor Ait is required to obtain the Gaussian cgs value 
from the implicit normalization of the magnetic field in 
the ZEUS code. When the surface strength of the dipole 
magnetic field is combined with Equation El it gives 

Ma,o 

the constant entropy, and m the number of baryons per particle, 
with inclusion of free electrons. For hot, completely ionized hydro- 
gen in the corona, m = 0.5, but in a cold disk m = 1. Si stands 
for the ideal gas constant 5R = 8.31 X 10 7 erg K — 1 mol - , from the 
ideal gas law p = pStT/m. 
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Fig. 3. — Initial angular velocity profiles along the disk mid-plane 
in our setups. Shown are the angular velocities in simulations with 
the stellar rotation rate Q* =0.3, 0.15 and 0.068, long-dashed line, 
dashed line and solid line, respectively. For comparison, the dotted 
line shows the Keplerian rotation profile. The kink in the initial 
rotational profile is near the initial inner disk radius Ri. 

where M 8 is the disk mass accretion rate in units of 
lCT 8 Af Q vr" 1 . For M 8 = 100 and M A o = 225, B* is 
about 100 G. 

2.2. Boundary conditions 

In order to mimic an absorbing stellar surface, we de- 
fine the "outflow" boundary condition aroun d a central 
object. As done bv lUchida fe Shibatal (j!985|) . we define 
part of the computational box surrounding the origin as 
a boundary layer, for a region further above the star. We 
set up a rotating circular layer of one grid cell thickness 
on top of the star, at a distance i?» from the origin, with 
the stellar rotation rate as a free parameter - see 
Figure [2j All the other hydrodynamical quantities are 
absorbed, so that the values in this layer are copied from 
the cells immediately above it. With this procedure, we 
neglect the stellar wind. 

The outer boundaries of our computational box are 
open, with the flows extrapolated beyond the boundary. 
One exception is a small part at the disk outer boundary. 
There, we prescribe a small mass inflow that is consis- 
tent with the initial radial component of the velocity in 
the disk. Reflection boundaries are imposed along the 
axis of symmetry and, in simulation SI, at the disk mid- 
plane inside the disk, where the normal component of the 
magnetic field is continuous, the tangential component is 
reflected, and the toroidal magnetic field is anti-reflected. 
Under axisymmetry for the disk mid-plane and the axis, 
B R (R = 0,Z) = B R {R, Z = 0) = 0, and with these con- 
ditions V • B — is satisfi ed. We use the Const rained 
Transport (CT) method of lEvans fc Hawlevl (1988) to en- 
sure that it is preserved to a machine round-off precision 
in computations. 

In simulation S2, we treat a small part of the disk mid- 
plane inside the disk gap as an open boundary, effectively 
extending the stellar absorbing layer into the disk gap. 
Such an extension is to ensure preservation of the disk 
gap in the simulation, even when the magnetic field is 
not strong enough to truncate the disk. This means that 



some other physical effects, like physical viscosity or ra- 
diation transfer, which we do not include in our simu- 
lations, would have to act in terminating the disk. Us- 
ing such a setup, we can study if a weaker stellar dipolc 
can launch outflows from the innermost magnetosphere. 
Caveat is that the final disk truncation radius in simu- 
lations S2 is then dependent on the initial setup, and is 
not self-consistently computed. 

2.3. Initial conditions 

We set up initial conditions for the density distribu- 
tion, velocity profiles, magnetic field and resistivity in 
the computational domain as follows. We set up a ro- 
tating disk that is simply gravitationally bound further 
away from the origin. For a self-consistent accretion disk, 
additional constraints such as constant fluxes through 
surfaces at different radii and stability to various modes 
of oscillation should be included. However, the disk sta- 
bility or the accretion process itself is not the subject of 
study here, and we treat the disk only as a supply of 
matter into the stellar magnetosphere. 

2.3.1. Density distribution 
The initial disk density distribution is 



p d (R,Z) 



R 



3/2 
off 



(R 2 oS + R 2 ) 3/i 



i{10~ 



1 - 



( 7 -l)Z 2 



2H 2 



1/(7-1) 



(11) 



shown in Figure [TJ The density is limited by the maxi- 
mum function, and ensured to be regular by a constant 
offset radius i? Q ff = 4. The disk is adiabatic with an in- 
dex 7 = 5/3, and physically thin, with an aspect ratio 
of H/R = 0.1, where H is the disk height at a given ra- 
dius R. For the initial inner disk radius we tried various 
initial positions of the initial inner disk radius Ri in our 
parameter study; here it is chosen to be at half size of 
the computational domain, Ri = 10i?* or half of that 
distance. It is not a critical parameter, as the disk will 
adjust its inner rim position during the simulation, but 
too close positioned Ri can, especially in a case of strong 
magnetic field, result in a too violent initial relaxation, 
which will stop a simulation. 

The corona above the star and in the disk gap coro- 
tates with the central object, and further away, with the 
underlying disk. The corona is in hydrostatic balance, 
with an initial coronal densitjU 



Pc(R,z) = 



(R 2 + Z 2 )- 3 / 4 



(12) 



which is obtained from the equality of gravitational and 
hydrostatic pressure. The free parameter k = Pd/Pc de- 
termines density in the corona. In similar studies n is 
usually in the range of 10 2 to 10 4 . In simulations with- 
out magnetic field and Sla we used k = 10 4 , and in Sib 
k = 10 . We address the influence of this parameter on 
the launching process in the resistive simulations, and 
investigate the range 10 2 to 10 6 for simulations S2. 

3 For such setup it is essential to set a f orce-free initial m a gnetic 
field in the computational box — see e.g. IFendt fc Elstncr iTl999l . 
l2000h . 
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2.3.2. Velocity profiles 

In our simulations, the initial stellar rotation rate is a 
free parameter, kept constant throughout the simulation. 
Since the time scale of change in stellar rotation is much 
longer than duration of simulations here, this constraint 
should not influence the outcome. In the case of T-Tauri 
type stars, there are observational i ndications that st ellar 
rotation rate is actually constant ()Irwin et al.ll2007j) . so 
that for those objects it is a plausible assumption even 
for very long lasting simulations. 

The corotation radius, at which matter in the disk is 
rotating with the angular velocity of the stellar surface, 
is: 



Rc 



= ( 



GAL 



1/3 



(13) 



The position of the corotation radius with respect to the 
disk truncation rad ius Rt defines two regi mes: R cor > 
R t , and R cor < R t . lUstvugova et al.1 (|2006t l and R09 in- 
vestigated the latter as a "fast rotating" (or "propeller" ) 
regime, and here we focus on the former, "slow rotating" 
regime. In this work we present results for a parame- 
ter study in a slow rotating regime, with stellar angular 
velocity of 0.15, which gives the rotation period of 11.8 
days. The corresponding corrotation radius is 10. li?*. 
For the disk, we adopt the following rotation profile: 



v^R, Z) = (1-e 2 



R 



1/2 
off 



{RIk + R 2 ) 1 '' 



■ exp 



The free parameter e gives the departure from the Ke- 
plerian rotation profile, and is chosen to be 0.1 in our 
typical simulations. For e = the disk would go back to 
the Keplerian profile. Figure [3] shows the initial angular 
velocity profiles at the equatorial plane of the disk. 

The initial poloidal velocity profile is given by the ra- 
dial inward velocity from accretion. The angular and 
sound speeds are both proportional to R~ x l 2 and, for a 
disk in hydrostatic equilibrium, the same is valid for the 
radial velocity. Both components of the poloidal velocity 
are given by 

v R (R, Z) = - ma e {R2 ^ R2y/i exp (-2^) , (15) 

v z (R,Z)=y R (R,Z)^ .(16) 

The constant parameter m s < 1 is used to obtain a sub- 
sonic inflow, and is chosen to be 0.1 in our simulations 
here. We also performed runs with m s = 0.3 and 0.6, 
which give larger influxes of mass into the disk, with 
similar results — for more massive disk, simulations are 
more prone to instabilities and tend to cease earlier than 
for lighter disk. The exponential factor in the equation 
effectively confines the initial disk profile. 

2.3.3. Magnetic field 

The initial magnetic field is a pure stellar dipole, and 
we computed it from the derivatives Br = —dA^/dZ 
and Bz = d(RA^,) /RdR of a magnetic potential: 



[i*R 



(R 2 + Z 2 fl 2 



(17) 
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Fig. 4. — The mass flux pv at T=300, in logarithmic color scale 
grading, in a simulation without magnetic field. There is no sig- 
nificant mass or angular momentum outflow after the relaxation 
phase. Vectors show the velocity of matter. The disk reaches the 
stellar surface, where matter is accreted onto the star. 

The stellar dipole magnetic moment /i* we set to unity. 
Setup with multipole expansion of magnetic field is feasi- 
ble in our simulations, but without stellar wind included 
we already neglect effects at the surface of the star. We 
assume that dipole is leading term in the disk gap and 
beyond. 

2.3.4. Resistivity and artificial viscosity 

The electrical resistivity r\ is defined through the elec- 
tric conductivity a as rj = c 2 /Aira, where c is the speed 
of light. The ratio of the advection and diffusion terms 
in the induction equation ([3]) is the magnetic Reynolds 
number, which is equal to the Lundquist number in our 
problem: 

Rm= lM^. (is) 
V 

The characteristic velocity is the Alfven speed «a,o = 
B /(47rp ) 1/2 at R (Z = 0). The reference time is the 
period of rotation at this radius Rq, ini>A,o> with «a,o ~ 

«K,0< 

To explain the physical processes, the accretion disk re- 
quires an enhanced, anomalous level of resistivity, which 
is much larger than the classical value. The anomalous 
resistivity could be an effect of MHD turbulence or am- 
bipolar diffusion in a partially ionized mediunQ. We 
set the initial constant resistivity of the disk to be of 
the same order of magnitude as the numerical resistivity 
rj = Ax 2 / At - 10~ 4 , which gives Rm ~ 10 4 , with At 
in units of rotation at the outer disk radius R ma x > where 
the time step is smallest. We find the numerical resistiv- 
ity by lowering the disk constant resistivity in the code 
until it does not affect the results. 

We omit the actual Ohmic term in the calculation of 
the MHD equations. When the Ohmic part is included in 
the internal energy equation in Zeus347, EquationUgains 
an additional Ohmic heating term —rjj 2 . The inclusion 
of this term is expensive computationally. However, the 
actual difference from the solution without the Ohmic 
term is negligible, as the pVv ter m is much large r . Sim - 
ilar results have been reported in IMiller k, Stone! ()1997f ) 
and R09. It would take the Ohmic term many orders of 
magnitude larger to produce a visible effect. 

4 For extensive di scussion of physical c onductivity in partiall y 
ionized disks see e.g. lWardle fc JNtJ 1 19991) , ISaTmeron et afl l2007h . 
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Fig. 5. — The mass flux pv for the case with medium magnetic 
field in the simulation Sla. Mass flux is shown in a logarithmic 
color grading, vectors show the velocity of matter, and the white 
solid line show the poloidal magnetic field lines. The disk reaches 
the stellar surface, where matter is accreted onto the star. A rar- 
efied, but fast axial outflow is launched from the magnetosphere 
immediately above the star. 

A resistive corona is essential for the magnetic re- 
connection to occur, and reconnection is crucial for 
re-organizing the magnetic field. The resistivity is 
modeled as a function of matter density, following 
iFendt fc Cernelul (pOOl . so that r) oc v A R ~ p^- 1 ^ 2 . 
For the adiabatic case 

V = r l0 p 1/3 ■ (19) 

To avoid unrealistically large r\ and too small Ohmic 
timesteps in the densest part of the domain, we limit 
this value to the order of unity, with t\q = 3.0. 

Another diffusive process in our simulations is the nu- 
merical viscosity. In our finite differencing numerical 
scheme, it is of the same order as numerical resistiv- 
ity. We do not treat the physical viscosity, only the 
von Neumann-Richtmyer artificial viscosity is included, 
through a constant parameter that controls the number 
of zones through which shocks are smoothed out. Such 
viscous term is significant only in presence of shocks. For 
a smooth flow, it is tiny, and for rarefactions, it is zero. 
The characteristic speed for viscous effects is the sound 
speed c s . 

3. SIMULATIONS 

We started with a disk in hydrostatic balance and, as 
a reference, performed a hydrodynamic simulation with- 
out magnetic field. The stellar rotation rate in the case 
shown in Figure |4] was set to f2* = 0.15, but we tried 
other f2*, smaller and larger, with the similar outcomes. 
After relaxation, the disk remained stable for hundreds 
of revolutions, in a quasi-stationary state, connected to 
the stellar equator, with the matter from the disk slowly 
falling onto the s tar. The disk got puffe d up, similar 
to the situation in iCemeliic fe FendD (j2004l ) , where there 
was no central object in the simulation, only the disk. 

3.1. Simulations with truncation of the disk by a strong 
magnetic field 

We seek to understand the effects of magnetic diffusion 
on the launching of outflows from the innermost vicinity 
of an YSO. In our code, numerical resistivity and numer- 
ical viscosity are of the similar order of magnitude. By 
including physical resistivity, but not physical viscosity, 
we probe the portion of parameter space with Pr m < 1. 
As mentioned in the Introduction, an YSO magnetic field 




Fig. 6. — Mass flux pv in logarithmic color grading for the case 
with strong magnetic field in the simulation Sib. The vectors show 
the velocity of matter, and the white solid lines show the poloidal 
magnetic field lines. Two outflows are launched from the magne- 
tosphere in the close vicinity of a star: a fast rarefied axial outflow 
immediately above the star, and a slower and denser conical outflow 
above the disk gap. The disk is truncated where its ram pressure 
and magnetic pressure are balanced. Poloidal magnetic field lines 
nearby the axis are removed from the shown sample, not to obscure 
the underlying mass flux distribution. 

is of the order of a few hundreds of Gauss to kG, so we set 
such magnetic field in our simulation. In the following 
text, we refer to such setup as simulations Sla. The stel- 
lar rotation rate remains the same as in the case without 
the magnetic field, £7* =0.15. 

The disk is falling radially onto the star, but with a dif- 
ference to the non-magnetic case that now a strong axial 
outflow is launched above the star. One example is shown 
in Figure The mass load in the outflow depends on 
the stellar rotation rate and the magnetic field strength, 
and on the disk accretion rate. In our axi-symmetric 
simulations it is not straightforward to conclude about 
the nature and stability of outflows launched co-axially, 
and so close to the central object defined as a bound- 
ary condition. We leave it for non-axisymmetric, full 3D 
simulations. 

To investigate realistic star-disk systems, we need a 
stable disk gap, which we do not obtain in any of simu- 
lations with magnetic field up to few hundreds of Gauss. 
Faster stellar rotation would help to establish it, because 
of a larger centrifugal force, but then the stellar rotation 
rate would become too large, compared to the range ex- 
pected for YSOs. We found that, with only the physical 
resistivity included, the only way to realistically obtain 
a disk gap is to increase the magnetic field. 

When we increase the magnetic field to the order of few 
hundred Gauss, magnetic pressure becomes sufficient to 
truncate the disk. When the disk gap stabilizes, in ad- 
dition to an axial outflow, a conical outflow is launched 
— see Figure [6l Such a result corroborates with simu- 
lations mentioned in fJTJ where, with an accretion rate 
of 10~ 8 M Q yr -1 , the magnetic field at which the disk 
becomes truncated is of the order of kG. We assign those 
simulations as Sib in the following text. 

Our purely resistive simulations are, therefore, repro- 
ducing the previously known results, that a sufficiently 
large magnetic field truncates the accretion disk — we 
discuss the truncation radius in greater detail in fj6] — 
and that a conical outflow is launched. With even 
stronger magnetic field and varying the accretion rate, 
it is possible to modify the gap extension and the inten- 
sity of both axial and conical outflows, but it is not clear 
if the reason for launching a conical outflow is a large 
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magnetic field, or other conditions near the disk gap. A 
caveat is that we neglected viscosity and radiation ef- 
fects, and hxed the stellar rotation rate. Therefore, it is 
not unexpected that, to truncate the disk, a large mag- 
netic field is needed. It remained the only adjustable 
parameter. 

In a setup as described, simulations with large mag- 
netic field tend to cease during, or not long after, the re- 
laxation, because of numerical problems. Without physi- 
cal viscosity, resistivity itself is not dissipative enough to 
stabilize the flow. To study the quasi-stationary state, we 
need simulations lasting for hundreds of rotations. Also, 
in order to anticipate including of the stellar outflow and 
radiative effects, which could help in establishing a gap, 
we need to devise a way of producing the disk gap with 
a smaller magnetic field, of the order of tens to hundreds 
of Gauss. 

3.2. Launching of a conical outflow 

To probe the portion of the parameter space with 
smaller stellar magnetic field, we devise a simulation in 
which the disk gap is numerically imposed, as described 
in i j2.2l Such simulation, dubbed S2 here, has been per- 
formed with a part of the disk mid-plane inside the disk 
gap defined as an open boundary. It means that the disk 
truncation radius is not determined self-consistentlj0- 

The results of our simulation S2 are shown in Figure[7J 
We show the density and poloidal mass flux pv p for the 
same time step in the right and left side of the same 
panel, to stress that in the density plots conical out- 
flows will typically not be visible even in logarithmic color 
grading. Instead, as shown in R09, outflows are well vis- 
ible in the poloidal mass flux plots in logarithmic color 
grading. This is probably one of the reasons why notice 
of the conical outflows was not made earlier, despite of 
many numerical efforts in the portion of parameter space 
where conical outflows should be appearing. Other fea- 
tures, as ejected plasmoid or accretion flow onto the star 
are well seen in both the density and mass flux plots, and 
have been described in the literature mentioned in Sjl] 

Now we obtained long lasting simulations, which in all 
respects resemble those from simulation Sib, but with 
the difference that they last longer, and the magnetic 
field required for launching of outflows is smaller for an 
order of magnitude. We can identify four evolutionary 
stages in progression in a system of an interacting mag- 
netosphere with its surrounding disk. In a case with 
a disk accretion rate 10 _7 M©/ year, for a rather small 
stellar dipole field of 38 G, the system goes through sim- 
ilar relaxation and initial evolution in all cases, with the 
similar geometry of the poloidal field. The results are 
robust in that they occur under a wide range of explored 
parameters, although each with different details. 

An initially pure dipole magnetosphere has already 
bulged out and brought some gas along with it at as 
few rotations as T — 2. Near the axis, some gas also 
flows out at high velocity due to magnetic pressure that 
is gradually building up, as shown in the large sizes of 

5 For a large enough magnetic field, of the order of kG as in simu- 
lation Sfb, such imposed disk gap is largely ignored by the disk, as 
matter is lifted above the disk equatorial plane. In the next paper, 
ICemeli ic & Shang (in preparation^, we describe configurations in 
which an accretion funnel onto the star can form. 




Fig. 7. — Snapshots in our simulation S2 with -R C or = 10-R*. 
To show the difference between density and mass flux plots, the 
left half of each panel shows the poloidal mass flux pv p , and 
the right half shows the density. Both plots are in logarithmic 
color scales, shown at the bottom of the panels. The poloidal 
magnetic field lines are shown in solid white lines, with lines 
in the axial region omitted, to show velocities in white arrows. 
The resolution is R x Z = (90 x 90) grid cells= (20 x 20)R„, 
and the initial density contrast between the disk and the corona 
is k = I0 4 . The initial magnetic field is a pure stellar mag- 
netic dipole, with B* = (380, 38, 3.8) G for a disk accretion rate 
(10 — 6 , 10 — 7 , 10 — 8 )M0/year, respectively. Each line shows, top to 
bottom, a characteristic stage discussed in this work: I) initial re- 
laxation when the magnetic field is swept in and pinched near the 
disk mid-plane toward the star, II) inflation and reconnection that 
end up opening the field, and strong infall of matter onto the star 
from the disk, III) retraction of the disk matter towards the coro- 
tation radius, with a transient inflow of matter onto the star, and 
the light bullets of fast matter expelled along the axis, IV) final, 
quasi-stationary stage, with a light, fast axial outflow and a dense, 
slow conical outflow. 
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Fig. 8. — A schematic sketch of the evolution in star-disk interaction. Stage one: the initial stellar dipole gets pinched during the 
relaxation, when matter flushes in toward the central star; Stage two: the magnetic field lines are open after reconnection, and disk matter 
can reach the surface of the star; Stage three: the disk matter retracts and a funnel flow forms from the disk inner radius and accretes a 
matter onto the star; Stage four: the system reaches a quasi-steady state, settling into a configuration consisting of an open stellar field 
and field footed in the disk. The arrows indicate the directions of the matter outflow. In the first three stages, the axial component is in 
some simulations strongly episodic (marked with gray shadow), on and off many times into the quasi-stationary state. In the fourth stage, 
a conical outflow forms, which is reaching a quasi-stationary state (marked with black shadowed arrow). 
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Fig. 9. — Change of the disk density profile in time, along the 
equatorial plane. We show densities at T=0, 5, 10, 100, 850 in thin 
solid, dotted, dashed, long dashed and thick solid line, respectively. 
The disk density is substantially modified only during the relax- 
ation, afterwards it does not change much. 
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Fig. 10. — Velocity in the axial direction in simulation S2, 
throughout the simulation. We show the velocity in the axial out- 
flow at (R, Z) = (OAR*, Z maI ) in a thin (red) solid line, and in the 
conical outflow at (R, Z) = (7.1-R* , Z max ) in a thick (black) solid 
line. 

the arrows there. The matter has flown in at a magnetic 
stagnation point around 8i?», where the magnetic field 
dragged in with the gas is pinched. Around T = 9, mat- 
ter went through a magnetic reconnection and ejected 
plasmoids. The reconnected and opened field enabled 
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Fig. 11. — Velocity profiles for results from the bottom panel in 
the Figure [7] at T=850. Profiles along the propagation direction 
of outflows in the Z-direction, parallel to the symmetry axis, at 
R = R t , we show in the ( Top panel) and at R = 5R* in the (Bot- 
tom panel). Components of the velocity in the Z, R and toroidal 
direction are shown in thin solid, dashed and dotted black lines, re- 
spectively, and the total velocity is plotted in the thick solid (black) 
line. Alfvcn velocity is plotted in the long-dashed (blue) line, es- 
cape velocity in (green) dot-short-dashed line, and the sound speed 
in dot-long-dashed (red) line. 

the disk gas to flush into the stellar surface and, at the 
same time, more violent gas flows are directed outwards 
both from the axial region and from the disk. At a later 
time, T = 45, after a few occurrences of the magnetic 
reconnection events, part of the field closes back to the 
stellar surface, and part remains open, footed near the 
new truncation radius. Matter channels through the field 
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Fig. 12. — Mass (Top panel) and angular momentum (Bottom 
panel) fluxes in the direction parallel to the axis of rotation, in a 
slice along the 7i m ax boundary, in the R-direction, in the quasi- 
stationary state from Figure [7] We show fluxes Fm and Fg at 
T=850 in solid line, and at T=852 in dashed line. The temporary 
influence of the disk is visible in the rightmost portion of the mass 
flux profile for T=850. It does not contribute to fast outflows, as 
the matter in the disk rolls back towards the disk. To avoid such 
apparent oscillation in the flux, we exclude fluxes from the outer 
portion of the box. Most of the angular momentum flux exits the 
system through the conical outflow, which also shows a periodic 
variation in intensity. 

lines that are footed both in the disk and the stellar sur- 
face. Matter along the axis is expelled in the form of 
bullets or, when more stabilized, in a more continuous 
stream. The bottom panel in Figure [7] shows the repre- 
sentative snapshot in those simulations at a much later 
time T = 850. The system settled into a configuration 
where the magnetic field has been opened into space with 
either foot in the star or in the disk, and formed loops 
that connect to both the star and the disk. The gas 
flows out from the axial regions on top of the star like a 
coronal wind and from the boundary on top of the loops 
along the diverging field lines open to the space, forming 
an outflow stream of conic shape. Although not driven 
completely by centrifugal mechanisms, the matter sur- 
rounding the axial region comes from the disk where one 
foot of the magnetic field is rooted. The disk material 
stays slightly outside of the magnetic footpoint where 
the field is pinched, around the truncation radius. When 
the disk accretion rate is well matched to the mass loss 
in the outflows, the simulations can last for hundreds of 
rotations. The disk, after the relaxation, does not differ 
much from the hydrodynamic case in Figure [4j 

What we observed in the simulations, and described 
in snapshots in Figure [3 taken at different times, seems 
to represent necessary steps in the evolution of the sys- 
tem consisting of a star with an accretion disk, when 
progressing towards a quasi-steady state. This process 
can be generalized into four conceptual stages, when the 
system attempts to evolve from an initial condition as a 



io-« - 




: 1 O- 3 — 



— 2x] O- 3 



Fig. 13. — Time evolution of mass (Top panel) and angular 
momentum (Bottom panel) fluxes in simulations S2, parallel to the 
axis of rotation, along the Zmax boundary. We show the fluxes 
after the relaxation in solid (black) line. In dashed (red) line we 
show the average value, computed starting from T=100, when the 
flow becomes quasi-stationary. The mass flux average value is 3.0 X 
10 -5 , and the angular momentum flux average value is 7.4 X 10 — 5 . 
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Fig. 14. — Average mass and angular momentum fluxes F m and 
Ft in our simulations with increasing density contrast ft between 
the disk and the corona. We show the mass flux in solid line, 
and the angular momentum flux in dashed line, in solutions with 
K = (10 2 , 10 3 , 10 4 , 10 5 , 10 6 ). Results for the mass flux for k = 10 2 
in our simulations are always larger for a factor of 2 to 3, than for 
larger, more realistic values of k. 

pure dipole threading everywhere into a disk in hydro- 
static balance. 

Stage I is the initial relaxation when the magnetic field 
is swept in and pinched near the disk mid-plane toward 
the star and the magnetic loops are twisted, inflating 
and forming plasmoids. The gas that flows with the 
field swirls in and is gradually accelerated in the axial re- 
gion by magnetic pressure being built up. Stage II takes 
the scene after the system goes through a reconnection 
that ends up opening the field, enabling strong infall of 
matter onto the star from the disk. The axial compo- 



Magnetospheric Accretion and Ejection of Matter in Resistive MHD Simulations 



11 




50 1 OO 1 
B. (Gauss) 



200 



FlG. 15. — Average mass and angular momentum fluxes F m and 
F( in our simulations with increasing stellar magnetic field. We 
show the mass flux in solid line, and the angular momentum flux in 
dashed line. With stellar magnetic field B* = (3, 10, 30, 77, 188)G, 
with disk accretion rate 10 — 6 MQ/year, both fluxes increase with 



increasing magnetic field. If we assume 10 
field strenghts are 1/10 of those values. 
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FlG. 16. — Change of resistivity with time in simulation S2. We 
show rj in two positions (R,Z) in the computational box: at the exit 
region of the outflow, in the position (7,Zmax), and in the middle 
of the computational box, in the position (10,10), in thick (black) 
and thin (red) solid lines, respectively. As expected in our model 
with resistivity in the disk corona dependent on the density, closer 
to the disk, resistivity is larger. 

nent strengthens, forming a series of bullets as a result 
of magnetic pressure from the twisted field. Stage III 
follows when the disk matter retracts towards the coro- 
tation radius, and a time-variable inflow of matter fun- 
nels onto the star from the inner disk truncation radius. 
The system may have several passages through the first 
three stages and finally move onto the quasi-steady state 
when the magnetic field is pinched and strong enough 
to balance the ram pressure of the disk gas, truncating 
the disk near the corotation radius. The magnetic field 
settles into a geometry where the field is open into the 
space both axially and conically, with some loops anchor- 
ing both in the star and the disk. Matter flows out along 
the open field lines, forming the axial and conical quasi- 
stationary outflows. Figure[8]is a schematic sketch of the 
stages. 

3.3. Properties of the conical outflow 

Here, we describe in more detail the properties of mat- 
ter in the quasi-stationary solutions with realistic mag- 



netic field dipole strength, in our simulations S2. 

Figure [9] shows the change of density in time along the 
disk equatorial plane. After the relaxation, our reservoir 
of matter in the simulation is not changing much. This 
means that the inflow of matter into the disk from the 
outer boundary, which mimics the accretion of matter 
from the outer part of the accretion disk, is well chosen. 
In the case of a too large inflow of matter into the disk, 
it would pile up and the disk would become unstable. 
On the other side, if the disk would become drained of 
matter, it would change the conditions we want to inves- 
tigate. 

In Figure [TU] we show the change of velocity in cho- 
sen points in the computational domain. Oscillations 
show a sign of instability working at a short time scale, 
so that both outflow components are not smooth even 
in the quasi-stationary phase. One possible reason is 
a permanent reshaping of the magnetic field by recon- 
nections along the flow — we discuss it in more detail 
in £j4j We show the velocity along the outflows in the 
quasi-stationary state in Figure 1111 The axial outflow 
is supersonic and sub-Alfvenic close to the stellar sur- 
face, but at the middle of the computational box it be- 
comes super- Alfvenic. Both the poloidal and total veloc- 
ity in conical outflows are larger than the escape velocity 
v esc = (2GM./R) 1 / 2 . 

We calculate the mass flux F m and the angular mo- 
mentum flux F^ in each half-plane, above or below the 
disk equator. They are defined as: 
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Figure Q2] shows the radial dependence of fluxes F m and 
Fi across the outer Z-boundary in the quasi-stationary 
state, and Figure Q2] shows the time evolution of those 
fluxes. To avoid influence of the disk in those computa- 
tions, we integrate fluxes only to R max /2 in each time 
step. During relaxation and stabilization of the system 
into quasi-stationary outflows, fluxes are large and vari- 
able, and become more steady afterwards. Reconnection 
events along the conical direction of the outflow con- 
tribute to the oscillation on the timescale of one rotation 
period. The contribution to the total fluxes from the 
axial component is two orders of magnitude smaller. 

3.4. Dependence on the density contrast and magnetic 
field strength 

Table [T] shows that simulations in the literature were 
performed for various contrasts of the disk to corona den- 
sity k. What is the influence of this parameter in our 
simulations? 

We show results of our parameter study in Figure UM 
There are no significant differences in the average fluxes 
in the range 10 3 to 10 6 . For k — 10 2 , the mass flux is 
always larger, so that it is double or, in some setups, 
triple the value from other cases. This trend, which is 
illustrated here for a case of sub-Keplerian rotation of 
the disk, is true also in the case of Keplerian rotation 
profile of the disk. Hence, results of simulations with 
k = 10 2 could be unrealistic in the case of launching of 
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Fig. 17. — Snapshots of poloidal magnetic field lines in the simulation S2 show the reshaping of the magnetic field during simulation. 
The initial stellar dipole is pinched into a plasmoid and ejected during relaxation process (Left panel). Afterwards, reconnection enables 
the change of magnetic field geometry (Middle panel) into the stellar and disk components (Right panel). 
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Fig. 18. — Elsasser number A in simulation S2. Above the line in 
the computational box where the Elsasser number A equals unity, 
it is A > 1 and launching of an outflow is possible. Below that line, 
A < 1 and the launching of matter can not occur. 

the astro-physical outflows. This is not surprising, taking 
into account that the realistic value of k is of the order 
of 10 s for YSOs. 

We also check how the mass fluxes in simulation S2 
depend on the magnetic field strength. Figure [T5] shows 
the average fluxes for magnetic fields up to the order 
of 0.2 kG, for YSOs with a disk mass accretion rate of 
10 _6 Mq yr -1 . Both the mass and angular momentum 
fluxes increase with increasing magnetic field strength. 

4. RECONNECTION AND OPENING OF MAGNETIC FIELD 

LINES 

Reconnection is essential for launching outflows from 
a star-disk magnetosphere. If, for some reason, recon- 
nection does not occur, a magnetic wall, or a magnetic 
tower form, and the evolution of the system will be differ- 
ent. The diffusive processes, as resistivity and viscosity, 
facilitate the reconnection, and shape the geometry of 
the final magnetic field. Even if physical resistivity is 
not included in the code, there is unavoidable numeri- 
cal resistivity, which can be estimated as 77 num = vAx, 
where Ax is the grid spacing. Numerical viscosity can be 
estimated the same way, with the coefficient ^ num - We 



included physical resistivity in the code. Its variation 
with time in two positions in the computational box is 
shown in Figure 1161 

In the simulation with non-resistive corona magnetic 
field does not relax smoothly, as it does in the simulation 
with a resistive corona. In a non-resistive simulation, the 
corona relaxes from the initial condition with a steeper 
gradient between the stellar and the disk component of 
the magnetic field, because of slower reconnection pro- 
cess. Applied to a stage II of our typical runs, it means 
that if a current sheet is persistent, as it would be in a 
high-resolution ideal-MHD simulations, a magnetic wall 
can form. The magnetosonic waves might reflect from it 
and destabilize the flow by strong shocks. 

In simulations with a disk gap, with non-negligible re- 
sistivity, the stellar dipole magnetic field is reshaped into 
the stellar and disk components. This occurs during the 
initial stages of the simulation and is later maintained or 
cyclically repeated. Figure [171 shows the reshaping of the 
magnetic field in simulations S2. It will look similar in 
any configuration of the star-disk setup. 

5. ELSASSER CRITERION FOR LAUNCHING RESISTIVE 
FLOWS 

Astrophysical outflows launched by a magnetospheric 
accretion-ejection are the result of the interaction of the 
active magnetosphere and the innermost portion of the 
disk. The Elsasser number can serve as a more general 
condition for launching resistive conical flows by mag- 
netic star-disk interactions. It is defined by the ratio of 
the Z component of the Alfven velocity and tj^Ik (see 
Sal meron et al.l (|2007t ) and references therein): 



A 



V 2 

V AZ 



> 1 



(21) 



In some systems, when the flow is mainly in the Z- 
direction, A can be similar or equal to the magnetic 
Reynolds number, but in general they are different. 
A > 1 is a better indicator of successful launching. In 
Figure 1181 we plot the Elsasser number A > 1 in simu- 
lation S2. We observe that A > 1 is valid only in the 
magnetosphere above the disk. For a small magnetic 
field of the order of 1 G, there is no region in the box 
that could satisfy A > 1, and no launching is possible. 
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Fig. 19. — A zoom into the closest vicinity of the star from 
Figure [6] in simulation SI, to show details in the disk gap, where 
the disk is truncated. Meaning of colors and lines is the same as 
in Figure [7] The stellar magnetic field is about 200 G, for the 
disk mass accretion rate 10 Mq yr -1 . Black arrow marks the 
position of the truncation radius of the disk, which is in this case 
Rt = 3.0R,. 

6. DISK TRUNCATION RADIUS 

The new truncation radius of the disk after the re- 
laxation is close to the line of balance of the disk ram 
pressure and the magnetic pressure, p + pv 2 — B 2 /8ir, 
located between a central object and the disk inner ra- 
dius — see Figure US In our results, it is close to the 
position where the disk density drops steeply, and where 
the magnetosonic Mach number equals unity at the equa- 
torial plane, as matter is being launched fast along the 
neutral line because of reconnection. The angular veloc- 
ity profiles along the disk mid-plane also show a dip at 
the disk truncation radii, passing from the gap into the 
disk. 

In simulations SI, different to setup in S2, we main- 
tained self-consistent boundary conditions at the disk 
mid-plane. Because of larger magnetic fields, such simu- 
lations lasted shorter than simulations S2, but enough 
that the disk would reach quasi-stationary state, and 
we can estimate the disk truncation radius from our 
results. In Figure [19] we show the disk truncation ra- 
dius in simulation Sib. It is measured by the position 
of the steep gradient in the disk density and it is at 
Rt = 3.0i?». The disk truncation radius can be esti- 
mated as R t = a t {BfRl 2 /2GM*M a 2 )i/7 , given by or- 
der of magnitude as the equilibrium of the ram pres- 
sure of a spherical envelope in free fall and the ma gnetic 
pressure of a stellar dipole (Eis ner fc Lambl l l977l ). The 
non-dime nsional factor at < 1 has bee n estimated to 
be 0.5 in IGhosh fc Lambl (fl978l H979aHbh and - 1 in 
IQstriker fc Shul (11995D . In our simulations the system 
departs significantly from the spherical infall case and 
results in a t for an order of magnitude too large, when 
compared to the estimate above. Simulation Sib yields 
a t = 5. 

iBessolaz et al.l (|2008l) estimate at to be equal to M s , 
where M s is the sonic Mach number measured at the 
disk mid-plane, using the radial velocity of the matter in 
the disk for a comparison with the sound speed. Their 
derivation was for the case with the accretion column 
onto the star, but conditions should be similar also in the 
case of outflow launched from the inner disk radius. For 
this estimate, inserting the values from simulation Sib, 



we obtain M s = 0.15, which gives a too small value of 
Rt = 0.3i?*. None of the two estimates is appropriate 
for our result here. 

A further investigation of location and stability of the 
disk truncation radius should be done in a simulation of 
the complete R-Z half-plane in cylindrical coordinates, 
to avoid effects of boundary conditions at the disk mid- 
plane. 

7. DISCUSSION 

In this work, we for the first time demonstrate launch- 
ing outflows with physical resistivity included in the 
whole computational box (and not only in the disk), 
without physical viscosity. We investigate effects of re- 
sistivity in the innermost vicinity of the central object 
surrounded by an accretion disk, and seek a common 
framework for similar works performed to date — some 
of them are shown in Table [TJ 

We set up our Zeus347 code without employing any 
special procedure to smooth the violent relaxation from 
initial conditions. With extensive exploration of param- 
eter space, the results, perhaps somewhat surprisingly, 
showed that star-disk systems in our purely resistive 
setup may undergo similar evolution to the results from 
the literature where additional free parameters or phys- 
ical processes, e.g. viscosity, were included. 

The initial evolution relies on two key ingredients: ini- 
tial closed magnetic loops that connect both the stellar 
surface and the innermost region of the surrounding disk, 
and sufficient resistivity in the system to enable the mag- 
netic reconnection to occur and for the field to continue 
its evolution to longer-lasting and more stable states. 

As sketched in Figure [SJ common phenomena have 
been identified in our simulations, which occur, some or 
all of them, in many simulations of interacting star-disk 
systems from Table[T] These phenomena can now be con- 
sidered robust after our extensive explorations in this pa- 
per. They are as follows. First, the system always relaxes 
as it transits from sometimes an unrealistic initial condi- 
tion to a more evolved configuration. The influx of mat- 
ter as the system relaxes pinches the magnetic field near 
the inner truncation radius. The twisting of the magnetic 
field lines that are connected to both the stellar surface 
and the inner parts of the disk inflates the loops. When 
enough twisting has been accumulated, magnetic recon- 
nection adjusts the topology of compressed field lines, 
and partially opens the originally closed loops. During 
this process, strong infall floods onto the stellar surface 
while the loops are opened up by the reconnection. After 
the reconnection, the loops reform and the field regains 
strength, with the flooded-in disk matter retracting back 
to the corotation radius. Irregular, transient funnels may 
occur as the system tries to adjust itself along the evo- 
lution to a longer-lasting state. In some setups, those 
processes can repeat themselves numerous times before 
a more or less quasi-equilibrium state is reached through 
the opening and closing of magnetic loops. This perhaps 
forms the basis of episodic accretion in the early phase 
of star-disk evolution. As for the later phase of an YSO 
evolution, there are also many cases in our simulations 
when the system switches into the final configuration im- 
mediately after the relaxation, and does not change any 
more. 

In this picture, the diffusion process, which is in our 
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case only the resistivity, plays an important role in ad- 
justing the system to a more sustainable configuration at 
hundreds of rotations. 

If there is a sufficiently strong stellar magnetic held, 
a two quasi-stationary outflow components form. Their 
intensity depends mainly on the rotation rates of the star 
and the disk, and on the dissipation processes. The mag- 
netic held lines are opened above the star, with the more 
or less episodic fast, light axial outflow, and the mag- 
netosphere between the star and the disk is the site of 
launching of the slower, conical component. Infall of mat- 
ter towards the axis might be obstru cted by a magnetic 
wall, if reconnection is suppressed (|Lvnden-Belll 119961 
2003), because resistivity can only moderately modify 
the flow shape on a slower time scale than the forma- 
tion of a magnetic wall. Such a wall could mount to a 
"magn etic tower" , reported in simulations bv lKato et al.1 
(2004) . More details on longer-lasting and stable funnels 
in our setup are being investigated in lCemeliic fc Shangl 
(|in preparation! ) . 

This axial, fast outflow component stabilizes in a last 
stage. High velocities along the axial region have often 
been obtained in simulations of astrophysical jets, but 
it was often considered a numerical artifact, and not a 
part of the magnetospheric accretion-ejection of matter 
in the star-disk system. H owever, it was present in m od- 
els based on observations ([Kwan fc Tademaru][l988l) and 
in models of evolution of an axial jet with reconnection of 
magnetic fi eld, related to disk evolution (|Goodson et al.l 
119971 11999T ). R09 also obtained a jet as a contribution 
to the outflow from the star-disk magnetosphere. We 
checked that our axial flow is dependent on conditions at 
the base of the flow (atop the star) and not formed just 
by magnetic collimation. The flow is highly variable at 
the beginning, as it builds up during the relaxation from 
the initial conditions, and during the intermittent phases 
of reconnection and inflation of the magnetic held, and 
remains variable in intensity throughout the simulation. 
A numerically produced flow would tend to be less depen- 
dent on actual physical parameters. However, to study 
the behavior of this outflow in greater detail, it would be 
needed to include stellar wind in the simulation, and to 
consider a larger physical domain in the computational 
box. Such a study would be more convincing, too, if per- 
formed in the full 3D simulations, as it could then answer 
questions about stability of the axial jet. We leave it for 
further study. 

In the first three transition stages, our results can 
be related to various stages obtained in many other 
simulations, with some of them l isted i n Table [H and 
presen ted also in lOstriker fc Shul (|1995l) , lLovelace et al.l 
(119951). The h i ghly ti me-dependent radial outflows in 
iGoodson et~aH ([ 19971 ); IGoodson fc Wingled (|1999f ) also 
resemble those we obtained here. Similar conclusions to 
ours can be found about the position of corotation ra- 
dius and to the occurrence of r econnection. For lo nger 
duration of simulations, which iKiiker et al.l (|2003l ) ob- 
tained for magnetic Prandtl number about unity, results 
about the disk truncation remain valid. The final stage 
in our simulations can be compared also with the results 
in R09, of which our simulations are necessarily a sub- 
set. Because of including physical viscosity, simulations 
of R09 have one degree of freedom more. A caveat is 



that their initial conditions differ from setups in other 
simulations in the literature, as well as ours, and the re- 
sults might be affected by this. We confirm that R09 
results are still valid, in the part concerning the resis- 
tive dissipation, with a more conventional initial setup. 
A quasi-stationary solution for the resistive simulation 
with substantial magnetic field in the star-disk magneto- 
sphere will consist of two outflow components, even when 
the viscosity is much smaller then resistivity. 

The strength of the stellar dipole, which we vary 
in the range of 0.1-200 G for the disk accretion rate 
lO _6 M0/year, affects the actual fluxes of mass and an- 
gular momentum. A very strong initial field may cause 
strong shear motion at the beginning, and can largely 
disturb the relaxation process, especially in the case 
of a non-resistive corona, because pinching and recon- 
nection of the magnetic field dep end critically on the 
conditions in the m agnetosphere ([Lovelace et al.l 119951 : 
IGoodson et al.l I1997T R In ideal MHD simulations per- 
forme d with a coarser grid — see e.g. iRomanova et abl 
(2002) — numerical resistivity can mimic resistive effects, 
and such results could be more realistic th an those in 
the h igh resolution ideal-MHD simulations ( Yan g et al.l 
1986). In high resolution simulations, physical resistivity 
should be included, or results could depend on numerical 
resistivity, which is not related to physical quantities of 
the setup. Any study of results is then necessarily flawed 
by purely numerical effects. However, a model for physi- 
cal res istivity itself ca n affec t the results, as was pointed 
out in IGoodson et~aH (|1997l ). 

In the accretion process, the angular momentum is 
partly transported through the system by means of hy- 
drodynamical viscosity. To some extent, viscosity is 
included in the effective resistivity when the magnetic 
Prandtl number is a ssume d to be of order unity — se e 
iFerreira fc Pelletierl (fl995h andlC asse fc Ferreiral (|2000j) . 
However, as shown in iMeliani et al.l (|2006f ). the viscous 
torque extracts angular momentum from the disk less 
efficiently than the magnetic torque for turbulent disks 
and can be neglected in simulations, as we have done in 
our simulations. 

To identify locations in which the launching is possi- 
ble in our computational box, we compute the Elsasser 
number A. For a sufficiently large magnetic field and 
resistivity, it shows that launching occurs from the mag- 
netosphere above the innermost disk region and the star. 
There is no disk component of outflow in our results. To 
self-consistently determine the contribution from the in- 
ner disk itself to the overall outflow from the star-disk 
system, larger portion of the disk and proper boundary 
conditions should be included in the simulation. Here we 
simulate only the innermost magnetosphere of the star- 
disk system, and outflow components which we obtain 
are then necessarily limited to the contribution from the 
disk gap. Essential role of this innermost part of the sys- 
tem is the re-shaping of the geometry of magnetic field 

6 It is worth noticing that the opening of the field lines does 
not happen only through reconnection. Another process which fa- 
cilitates it and can act in combination with reconnection is the 
inflation of the magnetic field lines. It occurs because differential 
rotation at the footpoints of the magnetic field loops which thread 
the star and t he disk, tends to ope n the field li nes. It has been de- 
scribe d in e.g. lGold & Hovlel plOT ) . lATvl l[T98fJ l and lLovelace et al.l 
(1995). 
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through the rcconnection, which is enabled by resistivity. 

In one of our setups, with large resistivity of the or- 
der of r\ = 0.5 in the disk, which is a thousand times 
larger than the numerical resistivity, matter is drained 
from the disk across the magnetic field lines. Such mech- 
anism could not lead to long-lasting outflows, except for 
an unusually large accretion rate, for example in close 
binaries. However, it is worth mentioning, as this could 
be one mechanism for disk disappearance in setups as 
ours. In astronomical objects, this could happen when 
the disk matter becomes stripped of some ingredients by 
planet formation or interaction with the environment, 
e.g. a neighbor star in close binary system. The resis- 
tivity in the vicinity of the central object could change, 
leading to an abrupt change in the accretion process, 
with the innermost part of the disk disrupted in time of 
the order of ten rotation periods at the inner disk ra- 
dius, much less than one rotation at the outer radius of 
the disk. Consequences could be observed in YSOs or at 
the end of the active accretion phase in planetary disks. 

8. SUMMARY 

We report results of our numerical simulations in the 
resistive MHD regime of magnetospheric accretion and 
ejection in the closest vicinity of a central object. For 
the first time, we show the launching of two long-lasting 
outflow components in the case of purely resistive mag- 
netosphere of a slowly rotating star. We also find that 
results, which are usually scattered in various setups and 
numerical methods in the literature, can be shown as four 
stages in the time-evolution of one setup for a star-disk 
system. This simplifies comparison of results from dif- 
ferent researchers. 

There are four common stages in our simulations: the 
flow relaxation with pinching of the magnetic field, the 
inflation with reconnection and partial opening of the 
stellar dipole field, the retraction of the disk towards the 
corrotation radius with transient funnel flow onto the 
stellar surface, and the equilibrium of magnetic and hy- 
drostatic pressure, with two outflow components from 
the innermost magnetosphere. The axial outflow is fast 



and of low-density, expelled by the magnetic pressure 
above the stellar surface. The conical outflow is slower 
and denser, and is launched by the magnetic force in the 
magnetosphere. 

In our simulations, only the stellar dipole is set as an 
initial condition for the magnetic field. The inclusion 
of a large scale interstellar or disk field does not change 
the outcome of simulations, as the stellar dipole is the 
leading term in the disk gap. The stellar dipole field 
also determines the disk inner radius, as it is positioned 
where the magnetic and the disk ram pressure balance. 
We identify possible location where the reconnection out- 
flow originates from the innermost magnetosphere by the 
Elsasser number A. The outflow components which we 
obtain in this work could just be part of the overall out- 
flow phenomena, or transient, when the inner portion 
of the disk actually participates dynamically and mag- 
netically. Current investigated solutions play a role in 
re-distribution of the initial stellar dipole magnetic field 
into the stellar and disk fields, enabled by resistivity. 

Our setup shares one caveat with most of the present 
works for a star-disk problem: we do not include the stel- 
lar wind in simulations. It would influence the solutions 
nearby the axis of symmetry and could even affect the 
very existence of the axial outflow component. Together 
with investigation of stability of axial outflows in full 3D 
treatment, we leave it for future work. 
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